function [V,M] = fastmvu(data, k, s)
%FASTMVU function [V,M] = fastmvu(data, k, s)
%       fast maximum variance unfolding
%
%       version 0.1 04/11/2014
%       Suqi Liu

disp('Loading data...');

X = double(data);
n = size(X,1);

disp('Performing KNN search...');

[IDX,D] = knnsearch(X,X,'K',k);

Ivec = reshape(repmat((1:n)',1,k-1),[],1);
IDXvec = reshape(IDX(:,2:k),[],1);
clear IDX;
Dvec = reshape(D(:,2:k).^2,[],1);
clear D;
Dvec = Dvec/max(Dvec);

A = sparse(Ivec,IDXvec,Dvec);
A = A + A';
A = (A~=0);

disp('Checking connectivity...');
%check connected
[S,~] = graphconncomp(A,'Directed',false);
if(S > 1) 
    error('Graph not connected! Try more neighbors.');
end

B = adj2inc(A);
clear A;

L = double(B)'*double(B);
clear B;
DD = diag(L);
L = 2*sparse(1:n,1:n,DD) - L;
clear DD;

disp('Finding eigenvalues...');

[V,~] = eigs(L,s+1,'SM');
clear L;
V = V(:,1:s);
Q = V(Ivec,:) - V(IDXvec,:);
clear Ivec;
clear IDXvec;

disp('Begin optimization...');

cvx_begin sdp
cvx_solver sedumi
    variable M(s,s) symmetric
    minimize (-trace(M))
    subject to
        sum(Q*M.*Q,2) <= Dvec
        M > 0
cvx_end
end